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Abstract We apply the theory of algebraic polynomials to analytically study the tran- 
sonic properties of general relativistic hydrodynamic axi symmetric accretion onto 
non-rotating astrophysical black holes. For such accretion phenomena, the conserved 
specific energy of the flow, which turns out to be one of the two first integrals of mo- 
tion in the system studied, can be expressed as a 8*'* degree polynomial of the critical 
point of the flow configuration. We then construct the corresponding Sturm's chain 
algorithm to calculate the number of real roots lying within the astrophysically rele- 
vant domain of M. This allows, for the first time in literature, to analytically find out 
the maximum number of physically acceptable solution an accretion flow with certain 
geometric configuration, space-time metric, and equation of state can have, and thus 
to investigate its multi-critical properties completely analytically, for accretion flow 
in which the location of the critical points can not be computed without taking re- 
course to the numerical scheme. This work can further be generalized to analytically 
calculate the maximal number of equilibrium points certain autonomous dynamical 
system can have in general. We also demonstrate how the transition from a mono- 
critical to multi-critical (or vice versa) flow configuration can be realized through 
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the saddle-centre bifurcation phenomena using certain techniques of the catastrophe 
theory. 

Keywords accretion, accretion discs • black hole physics • hydrodynamics • 
gravitation 

1 Introduction 

In order to satisfy the inner boundary conditions imposed by the event horizon, ac- 
cretion onto astrophysical black holes exhibit transonic properties in general [1]. A 
physical transonic accretion solution can mathematically be realized as critical so- 
lution on the phase portrait [2, 3, 4,5, 6, 7 ,8 ,9, 10, 1 1, 12, 13|. Multi-critical accretion 
may be referred to the specific category of accretion flow configuration having multi- 
ple critical points accessible to the accretion solution. For certain astrophysically rele- 
vant values of the initial boundary conditions, low angular momentum sub-Keplerian 
axisymmetric black hole accretion can have at most three critical points all together 
- where two saddle type critical points accommodate one centre type critical point 
in between them [1 , 14, 15, 16, 17, 18, 19,20,21 ,22,23,24,25,26,27,28,29,30,31 ,3j 
[33l[34l . Transonic solution passing through the aforementioned two critical points can 
be joined through a stationary shock generated as a consequence of the presence of 
the angular momentum barrier [IJ. 35. 29. 36. 32. 33]. The existence of such weakly 
rotating accretion in realistic astrophysical environment have also been observed ll37l 
[38l[39ll40ll4ni42 ! . A complete investigation of the multi-critical shocked accretion 
flow around astrophysical black holes necessitates the numerical integration of the 
nonlinear stationary equations describing the velocity phase space behaviour of the 
flow. 

However, for all the importance of transonic flows, there exists as yet no gen- 
eral mathematical prescription allowing one a direct analytical understanding of the 
nature of the multi-criticality without having to take recourse to the existing semi- 
analytic approach of numerically finding out the total number of physically accept- 
able critical points the accretion flow can have. 

This is precisely the main achievement of our work presented in this paper Using 
the theory of algebraic polynomials, we developed a mathematical algorithm capable 
of finding the number of physically acceptable solution a polynomial can have, for 
any arbitrary large value of n (n being the degree of the polynomial). For a specified 
set of values of the initial boundary conditions, we mathematically predict whether 
the flow will be multi-critical (more than one real physical roots for the polynomial) 
or not. This paper, thus, purports to address that particular issue of investigating the 
transonicity of a general relativistic flow structure around non rotating black holes 
without encountering the usual semi-analytic numerical techniques, and to derive 
some predictive insights about the qualitative character of the flow, and in relation 
to that, certain physical features of the multi-criticality of the flow will also be ad- 
dressed. In our work, we would like to develop a complete analytical formalism to 
investigate the critical behaviour of the general relativistic low angular momentum 
inviscid axisymmetric advective hydrodynamic accretion flow around a non rotating 
black hole. 
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To accomplish the aforementioned task, we first construct the equation describing 
the space gradient of the dynamical flow velocity of accreting matter. Such equation 
is isomorphic to a first order autonomous dynamical system. Application of the fixed 
point analysis enables to construct an 8th degree algebraic equation for the space 
variable along which the flow streamlines are defined to possess certain first integrals 
of motion. The constant coefficients for each term in that equation are functions of 
astrophysically relevant initial boundary conditions. Such initial boundary conditions 
span over a certain domain on the real line K - effectively, as individual sub-domain 
of RxMxR for the polytropic accretion. The solution of aforesaid equation would 
then provide the critical (and consequently, the sonic) point rc- The critical points 
itself are permissible only within a certain open interval ]rg, L^oo[, where rg is the 
radius of the event horizon and L_>oo is the physically acceptable maximally allowed 
limit on the value of a critical point. 

Since for polynomials of degree n > A, analytical solutions are not available, 
we use the Sturm's theorem (a corollary of the Sylvester's theorem), to construct the 
Sturm's chain algorithm, which can be used to calculate the number of real roots 
(lying within a certain sub-domain of M) for a polynomial of any countably finite 
arbitrarily large integral n, subjected to certain sub-domains of constant co-efficients. 
The problem now reduces to identify the polynomials in rc with the Sturm's se- 
quence, and to find out the maximum number of physically acceptable solution an 
accretion flow with certain geometric configuration, space-time metric, and equa- 
tion of state can have, and thus to investigate its multi-critical properties completely 
analytically, for accretion flow in which the critical points can not be computed an- 
alytically. Our work, as we believe, has significant importance, because for the first 
time in the literature, we provide a purely analytical method, by applying certain the- 
orem of algebraic polynomials to check whether certain astrophysical hydrodynamic 
accretion may undergo more than one sonic transitions. 

We further demonstrate how the transition of number of critical points may be 
taken into account considering the bifurcation phenomenon in the parameter space. 
The transition of number of critical points in this case is associated with the merging 
and destruction (or emergence and separating apart, viewing in the other way round) 
of a saddle-centre pair, i.e. a saddle-centre bifurcation common in conservative sys- 
tems, which may be tracked down using technique of catastrophe theory. The bifur- 
cation lines in the parameter space exactly conform with the transition boundaries of 
the across which the number of critical points changes. 



2 First Integral of Motion as a Polynomial in Critical Radius 

Following standard literature, we assume that the axisymmetric accretion flow has a 
radius dependent local thickness H{r), and its central plane coincides with the equa- 
torial plane of the black hole. It is common practice in accretion disc theory (|43 ,44] 
l451l46ll26ll47ll48l|49ll50ll5Ti ) to use the vertically integrated model in describing the 
black hole accretion discs where the equations of motion apply to the equatorial plane 
of the black hole assuming the flow to be in hydrostatic equilibrium along transverse 
direction. We follow the same procedure here. The thermodynamic flow variables are 
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averaged over the disc height, i.e. a thermodynamic quantity y used in our model is 
vertically integrated over the disc height and averaged as y = /o^''^'' ydh/ jl^^^ ' dh. 

We follow 1,5 2 J to derive an expression for the disc height H{r) in our geometry 
since the relevant equations in 11521 are non-singular on the horizon and can accom- 
modate both axial and quasi spherical flow geometry. The disc height comes out to 
be I32J, 



/2(7-l)(l-M2)[r3-A^(r-2)] 
7[7-(l + c2)](r-2) 



(1) 



where A and 7 are the specific flow angular momentum and the adiabatic index of 
the flow, respectively, u and Cg being the dynamical flow velocity and the speed of 
propagation of the acoustic perturbation (adiabatic sound speed) embedded within 
the accretion flow. In this work, we employ polytropic accretion. However, poly- 
tropic accretion is not the only choice to describe the general relativistic axisymmet- 
ric black-hole accretion. Equations of state other than the adiabatic one, such as the 
isothermal equation |27 1 or two temperature plasma ll53l have also been used to study 
the black-hole accretion flow. 

For accretion flow of aforementioned category, two first integrals of motion along 
the streamline, viz, the dimensionless conserved specific flow energy i.e., the energy 
per unit mass which actually is scaled by the rest mass of the flow 8, and the mass 
accretion rate M, may be obtained as (the radial distance r here is actually scaled by 
the factor GMbh /c^, and all the velocities, both u as well as Cg have been scaled 
by the velocity of light c in vacuum. Mbh is the mass of the black hole. Natural 
geometric unit has been used where the values of all fundamental constants have 
been taken to be unity, see, e.g., |[32l for further detail) 



(7-1) 
7-(l + c2) 



r-2 



1 



(2) 



M 



inpcsr^u /2(7-l)[r3-A2(r-2)] 
A V 7[7-(l + c2)] 



(3) 



where p is the mass density. The expression for £ is obtained by integrating the sta- 
tionary part of the Euler equation and the expression for M is obtained by integrating 
the stationary part of the continuity equation (by properly taking care of the flow 
thickness). The conserved specific entropy accretion rate ^A is computed as a quasi 
constant multiple of M as: 



M 



1 



7-1 



AV(r 



(4) 



We thus have two primary first integrals of motion along the streamline - the specific 
energy of the flow £ and the mass accretion rate M. Even in the absence of creation 
or annihilation of matter, the entropy accretion rate A4 is not a generic first integral 
of motion. As the expression for Ai contains the quantity K=p/ p'' (p being the flow 
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pressure), which is a measure of the specific entropy of the flow, the entropy accretion 
rate Ai remains constant throughout the flow only if the entropy per particle remains 
locally invariant. This condition may be violated if the accretion is accompanied by 
a shock. Thus A4 is conserved for shock free polytropic accretion and becomes dis- 
continuous (actually, increases) at the shock location, if such a shock is formed. 
The gradient of the acoustic velocity Cg as well as the dynamical velocity u can be 
obtained by differentiating the expression for the entropy accretion rate and the mass 
accretion rate respectively: 



dr 



,(7-1) b-{l + cl)] 



where 



(7- 



/i(r,A) 



u dr 



3r^ 



r4 - A2r(r - 2) 



du 
dr 



{^)clh{r. A) - /2(r, A) ^ AA(r,A,c.) 

I?(u,Cs) 



2cl 
"(7+1) 



where 



/2(r-,A) = 



2r - 3 



(5) 



(6) 



(7) 



(8) 



r(r - 2) 7-4 - A2r(r - 2) ' 

A real physical transonic flow must be smooth everywhere, except possibly at a 
shock. Hence, if the denominator 2?(u,Cs) of Eq. O vanishes at a point, the nu- 
merator A/^ (r, A, Cs) must also vanish at that point to ensure the physical continuity 
of the flow. One therefore arrives at the critical point conditions by making T) {u, Cg) 
and (r, A, Cg) of Eq. ^ simultaneously equal to zero. We thus obtain the critical 
point conditions as 



h{r,,\) + f2{r,,\)' 



h + 1 



/2(rc,A) 
/i(r-c,A) 



(9) 



where Uc = u{rc) and Cc = Cs{rc), r^ being the location of the critical point. 
/i(rc, A) and /2(7'c, A) are defined as: 



/i(rc,A) 



3r; 



2A2rc + 3A2 



A2r,(r,-2) 



/2(re,A) 



2rr 



2rl 



AV, + A2 



rc(rc-2) r^-X^r^ir 



2) 
(10) 

Clearly, the critical points are not coincident with the sonic points since Mc — 
(uc/cc) < 1. This is a consequence of the choice of the equation of state. The adi- 
abatic equation of state used in this work produces non constant (with respect to the 
radial space direction) sound speed. Since the disc height contains the sound speed 
and the thermodynamic quantities calculated in the accretion flow have been averaged 
over the flow thickness, non constant sound speed accounts for the non-isomorphism 
of the critical points and the sonic points. If one uses the sound speed obtained from 
isothermal equation of state, or a flow geometry different from the configuration in 
the vertical equilibrium as has been assumed here, the critical points will coincide 
with the sonic points, see, e. g., 1131113411 for further detail. 
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We substitute the explicit value of Uc and Cc from Eq. (|9]i to the expression for 
the specific energy £ in Eq. (|2]i to derive the explicit form of the energy first integral 
polynomial in Tc as: 

r^{-36 (-1 + (-1 + 7)^} + rj{12 (-1 + 7) (-17 (-1 + 7) + (-H + 137))} + 
r^{-24(-l + 7)^ (-16 + A^) + f 2 (-121 + GOA^ + 7 (286 - m\^) + 7^ (-169 + SGA^))} 
+rl{-2 (120 + (-86 + 163£2) + 7^ (120 + (-86 + 99£2) A^) - 27 (120 + (-86 + liiS'^) A^))} 

(^-460(-l + -if + f 2 (588 - 2b\^ + 7^ (356 - gA^) + 7 (-976 + SOA^)))} 

+r^{4A2 (^136(-1 + 7)^ + £2 (-88 + 45A2 + 7 (148 - 52A2) + 7^ (-52 + ISA^))) } 

+r2{-4A2 (60 + 121£2a2 + ^2 (go + ilE'^X^) - 27 (60 + QTE^X^))} 

+rc{32£2 (18- 197 + 57^) A*} + {-64f2(-2 + 7)^A"*} = (11) 

The above equation, being an n = 8 polynomial, is non analytically solvable. Being 
equipped with the details of the Sturm theorem and its appropriate application in the 
next section (^J3]l, in f|4] we will demonstrate how we can analytically find out the 
number of physically admissible real roots for this polynomial, and can investigate 
the transonicity of the flow. 

3 Sturm theorem and generalized sturm sequence (chain) 

In this section we will elaborate the idea of the generalized Strum sequence/chain, 
and will discuss its application in finding the number of roots of a algebraic poly- 
nomial equations with real co-efficients. Since the central concept of this theorem is 
heavily based on the idea of the greatest common divisor of a polynomial and related 
Euclidean algorithm, we start our discussion by clarifying such concept in somewhat 
great detail for the convenience of the reader. 



3.1 Greatest common divisor for two numbers 

Given two non-zero integers zi and Z2, one defines that zi divides Z2, if and only 
if there exists some integer zaG Z such that: 

Z2 = 23^1 (12) 

The standard notation for the divisibility is as follows: 

zi I Z2 rneans 'zi divides 22' (13) 

The concept of divisibility applies to the polynomials as well, we treat such situations 
in the subsequent paragraphs. 

Now consider two given integers zi and Z2, with at least one of them being a non- 
zero number The 'greatest common divisor' (or the 'greatest common factor' or the 
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'highest common factor') of zi and Z2, denoted by gcd{zi, is the positive integer 
Zd'E.'i, which satisfies: 

\)zd\zi and Zd|^2- 
ii)For any other z^^'L^ if Zc\zi andzc|^2 

thenzcl^d (14) 

In other words, the greatest common divisor gcd{zi, Z2) of two non zero integers 
Z\ and Z2 is the largest possible integer that divides both the integers without leaving 
any remainder Two numbers zi and Z2 are called 'co-prime] (alternatively, 'relatively 
prime'), if: 

gcd{zi,Z2) = l (15) 

The idea of a greatest common divisor can be generaUzed by defining the greater 
common divisor of a non empty set of integers. If Sz is a non-empty set of integers, 
then the greatest common divisor of Sz is a positive integer Zd such that: 

i) If Z(j|zifor all zi^Sz 
ii)If Z2|zi, for all ziG^z, thenz2|Z(i (16) 

then we denote Zd = gcd(Sz)- 



3.2 Euclidean algorithm 

Euclidean algorithm (first described in detail in Euclid's 'Elements' in 300 BC, and 
is still in use, making it the oldest available numerical algorithm still in common 
use) provides an efficient procedure for computing the greatest common divisor of 
two integers. Following Stark f54l, below we provide a simplified illustration of the 
Euclidean algorithm for two integers: 

Let us first set a 'counter' i for counting the steps of the algorithm, with initial 
step corresponding to z = 0. Let any zth step of the algorithm begins with two non- 
negative remainders ri_i and ri_2 with the requirement that r^-i < ri_2, owing to 
the fact that the fundamental aim of the algorithm is to reduce the remainder in suc- 
cessive steps, to finally bring it down to the zero in the ultimate step which terminates 
the algorithm. Hence, for the dummy index i, at the first step we have: 

r_2 = Z2 and r_i = zi (17) 

the integers for which the greatest common divisor is sought for After we divide Z2 
by Zi (operation corresponds to i — 1), since Z2 is not divisible by zi, one obtains: 

r-2^qor-i+ro (18) 

where tq is the remainder and qq be the quotient. 

For any arbitrary ith step of the algorithm, the aim is to find a quotient qj and 
remainder r^, such that: 



n^2 = qin-i + n, where ri < ri_i 



(19) 
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at some step i — j (common sense dictates that j can not be infinitely large), the 
algorithm terminates because the remainder becomes zero. Hence the final non-zero 
remainder 7 j_i will be the greatest common divisor of the corresponding integers. 

We will now illustrate the Euclidean algorithm for finding the greatest common 
divisor for two polynomials. 

3.3 Greatest common divisor and related Euclidean algorithm for polynomials 

Let us first define a polynomial to be 'monic' if the co-efficient of the term for the 
highest degree variable in the polynomial is unity (one). Let us now consider pi (x) 
and p2 {x) to be two nonzero polynomials with co-efficient from a field F (field of real, 
complex, or rational numbers, for example). A greatest common divisor of pi {x) and 
P2 {x) is defined to the the monic polynomial pd [x) of highest degree such that pd [x) 
divides both pi{x) and P2{x). It is obvious that F be field and Pd{x) be a monic, are 
necessary hypothesis. 

In more compact form, a greatest common divisor of two polynomials pi , p2GK[X 
is a polynomial p^gRM of greatest possible degree which divides both pi and p2- 
Clearly, pd is not unique, and is only defined upto multiplication by a non zero scalar, 
since for a non zero scalar cgK, if pd is a gcd{pi,p2'^^\^), so as cpd- Given poly- 
nomials pi,p2SK[X], the division algorithm provides polynomials p3,p4SR[X], with 
deg{p4) < deg{p3) such that 

Pi = P3P2 + Pi (20) 

Then, if pd is gcd{pi,p2), if Sind only if pd is gcd{p2,P4:) as is obvious. 

One can compute the gcd of two polynomials by collecting the common factors 
by factorizing the polynomials. However, this technique, although intuitively simple, 
almost always create a serious practical threat while making attempt to factorize the 
large high degree polynomials in reality. Euclidean algorithm appears to be relatively 
less complicated and a faster method for all practical purposes. Just like the integers 
as shown in the previous subsection, Euclid's algorithm can directly be applied for 
the polynomials as well, with decreasing degree for the polynomials at each step. The 
last non-zero remainder, after made monic if necessary, comes out to be the greatest 
common divisor of the two polynomials under consideration. 

Being equipped with the concept of the divisibility, gcd and the Euclidean al- 
gorithm, we are now in a position to define the Strum theorem and to discuss its 
applications. 

3.4 The Sturm Theorem: The purpose and the definition 

The Sturm theorem is due to Jacaues Charles Francois Sturm, a Geneva born French 
mathematician and a close collaborator of Joseph Liouville. The Sturm theorem, pub- 
lished in 1829 in the eleventh volume of the 'Bulletin des Sciences de Ferussac' under 
the title 'Memoire sur la resolution des equations numeriques'0. The Sturm theorem. 



' According to some historian, the theorem was originally discovered by Jean Baptist Fourier, well 
before Sturni, on the eve of the French revolution. 
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which is actually a root counting theorem, is used to find the number of real roots over 
a certain interval of a algebraic polynomial with real co-efficient. It can be stated as: 

Theorem 1 The number of real roots of an algebraic polynomial with real coefficient 
whose roots are simple over an interval, the endpoints of which are not roots, is equal 
to the difference between the number of sign changes of the Sturm chains formed for 
the interval ends. 

Hence, given a polynomial pGK[X], if we need to find the number of roots it can have 
in a certain open interval ]a, b[, a and b not being the roots of /, we then construct a 
sequence, called 'Sturm chain', of polynomials, called the generalized strum chains. 
Such a sequence is derived from p using the Euclidean algorithm. For the polynomial 
p as described above, the Sturm chain po,pi... can be defined as: 

Po=P 
Pi ^P' 

p„ = -rem(p„_2,p„-i) ,"->2 (21) 

where rem {pn-2,Pn-i) is the remainder of the polynomial p„_2 upon division by 
the polynomial Pn-i- The sequence terminates once one of the pi becomes zero. 
We then evaluate this chain of polynomials at the end points a and b of the open 
interval. The number of roots of p in ]a, b[ is the difference between the number of 
sign changes on the chain of polynomials at the end point a and the number of sign 
changes at the end point b. Thus, for any number t, if Np^^f^ denotes the number of 
sign changes in the Sturm chain po{t),pi{t), then for real numbers a and b that 
(both) are not roots of p, the number of distinct real roots of p in the open interval 
]a,b[ is [Npi^aj ~ ^p{b)] - By making a->~oo and 6— >+oo, one can find the total 
number of roots p can have on the entire domain of K. 

A more formal definition of the Strum theorem, as a corollary of the Sylvester's 
theorem, is what follows: 

Definition Let R be the real closed field, and let p and P be in i?[X].The Sturm 
sequence of p and P is the sequence of polynomials {po,pi, ■■■,Pk) defined as fol- 
lows: 

Po =p,Pi^ p'P 

Pi = Pi-iqi - Pi-2 with qi E R[X] and deg(pi) < deg{pi-i) for i = 2, 3, k, 
Pk is a greatest common divisor of p and p'P. 

Given a sequence (ao, ak) of elements of R with qq ^ 0, we define the number 
of sign changes in the sequence (ao, ...ak) as follows: count one sign change if 0^0; < 
with l = i + loTl>i + l and aj — for every j, i < j < I. 

If a € i? is not a root of p and {po, ...,pk) is the Sturm sequence of p and P, we 
define v{p, P; a) to be the number of sign changes in {po{a), ...pk{a)). 

Theorem 2 (Sylvester's Theoren^l) Let Rbe a real closed field and let p and P be 

two polynomials in R[X]. Let a,b E R be such that a < b and neither a nor b are 
roots of p. Then the difference between the number of roots of p in the interval ]a, b[ 
for which P is positive and the number of roots of p in the interval ]a, b[for which P 
is negative, is equal to v{p, P; a) — v{p, P\ b) 



^ As stated in (55]. 
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Corollary 1 (Sturm's Theorem): Let R be a real closed field and p G R[X]. Let 
a,b ^ R be such that a < b and neither a nor b are roots of p. Then the number of 
roots ofp in the interval ]a, b[ is equal to v(j), 1; a) ~ v{p, 1; b). 

The proof of these two theorems are given in the Appendix I. 

4 Number of available critical points for relativistic accretion 

We first write down the complete expression for the Sturm chains. Then for a suitable 
parameter set [£, A, 7], we can find the difference of the sign change of the Sturm 
chains at the open interval left boundary, i.e., at the event horizon and at the right 
boundary, i.e., at some suitably chosen large distance, say, 10^ gravitational radius 
(which is such a large distance that beyond which practically no critical point is ex- 
pected to form unless the specific flow energy has an extremely low value, i.e., very 
cold accretion flow), to find the number of critical points the accretion flow can have. 

The form of the original polynomial has already been explicitly expressed using 
left hand side of Eq.[TT] We now construct the Sturm chains as: 

Po{r) = asr^ + a-jr"^ + a^r^ + asr^ + a4r'^ + a^r^ + a2r^ + air + oq 
pi{r) = Sagr'^ + lajr^ + Qa&r^ + ba^r^ + 4a3r^ + Sasr^ + 2a2r + ai 
P2{r) = -rem{po/pi) = cer^ + c^r^ + c^r'^ + csr^ + C2r^ + cir + cq 
(the negative of the remainder of division of po by pi) 



P3{r) = 


-rem{pi/p2) 


- d5r^ 


+ d4r* 


+ d3r^ 


+ d2r^ + dir + do 


P4{r) = 


- rem (P2/P3) 


4 

= - 




-f e2r^ - 


f eir + eo 


P5{r) = 


~~rem{p3/p4) 


- fsr^ - 




+ hr + 


/o 


Peir) = 


-rem{pi/pz) 


= 92r^ ^ 


f 9ir 4 


90 




pr{r) = 


-rem{p5/p(i) 


= hir^ 


- K 






Psir) = 


-rem{pe/p7) 











Where the explicit expression of the corresponding co-efficients ai,Ci,di... has been 
provided in the equation (fTTl l and in the Appendix - II. If one needs to figure out 
the number of roots of po in [a,b], the number of sign changes in the sequence 
Po{a),pi{a),p2{a),p3{a),p4,{a),p5{a),pe{a),pT{a),ps{a) is to be counted and let 
us call it v{po, a). Similarly, the count the number of sign changes in the sequence 
Po{b),Pi{b),P2{b),P3{b),p4{b),P5{b),pe{b),p7{b),ps{b) is to be called as v{po,b). 
Then, the number of roots of po in [a, b] is v{pq, a) — v{po, b). 

It is important to note that direct application of the Sturm's theorem may not al- 
ways be sufficient since some of the roots may yield a negative energy for £ (since the 
£ equation was squared to get the polynomial). Since we are interested in accretion 
with the positive positive Bernoulli's constant, to get positive values of the energy, 
we must impose the condition that 



7-(l + cD >0, 



(22) 
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which is the term present in £ which could go negative. This introduces the condition 

that ——r- > 0, where p(r) and q(r) are 4th order polynomials given by, 

q[r) 

p{r) = 6(7 - - (II7 - 13)r^ 

-(57 - 3)X'^r^ + 2(97 - 5)X'^r - 8(27 - IjA^, (23a) 
q(r) = 6r^ ~ 12r^ - AX^r^ + 14AV - 12X^. (23b) 

To find the region where this happens, one has to find the 4 roots of each of p{r) and 
q{r) - which is analytically possible since roots of quartics are analytically solvable. 
Once the roots are obtained it is a trivial matter to check for what regions the rational 
function is positive. 




Fig. 1 The lighter region (online version red) corresponds to 3 roots and the darker shade (online version 
blue) indicates 1 root only. The value of 7 is 4/3. 



A simplified version for the above mentioned procedure to find the positivity 
condition is as follows: 

We would like to find out the intervals in which p{r)/q{r) > where p{r) and 
q{r) are quartic polynomials. We factorize p{r) = {r — ri){r — r2){r — r3){r — r^) 
and q{r) = (r — si)(r — S2)(r — S3)(r — S4) using the algorithm for finding roots 
of a quartic. If the roots are all real, we note down the sign changes of each factor to 
the right and left of each root and find out the intervals where the rational function 
is positive. If there are complex roots, they come in complex conjugates, since the 
coefficients of the polynomials are real. Say, if is complex and r4 is its complex 
conjugate, then the part (r — r'3)(r — r4) = — {r^ + r4)r + r^r^ does not change 
sign since it is non-zero on the real line. It is easy to determine its sign. 

To demonstrate the procedure described above, the number of roots of the 8th 
order polynomial po (in the Strum sequence) within the admissible range of f , A and 
7 (usually by keeping the value of 7 to be fixed to obtain a two dimensional parameter 
space) are evaluated explicitly and that shows two distinct regions in £ — A space (see 
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Fig ID- The wedge shaped region corresponds to 3 roots implying 3 critical points and 
the rest of the parametric space corresponds to single root implying only one critical 
point. This feature emerging from the above mentioned algorithm, exactly conforms 
with the numerical results (using the explicit root finding methods) available in the 
current literature [32 1. It may be worthwhile to mention here that in addition to these 
roots there exists another root for the the whole range of parameter space shown 
in the Fig. [T] that is located very near to the event horizon (i.e. within 1-1.5 times 
Schwarzchild radius), but being a centre it is physically untenable to be a sonic point 
(a critical point through which a physical accretion solution, connecting the event 
horizon with to infinity, can pass) and hence has always been justifiably ignored in 
the literature. 



1.008 
1.006 

8 

1.004 
1.002 
1 

2.6 2.8 3 3.2 3.4 3.6 3.8 4 

X 

Fig. 2 Boundary of transition: Contour line det (3) = (for 7 = 4/3). 

The transition boundaries from ni number of roots to 712 number of roots, in the 
parameter space, can be more easily obtained using catastrophe theory. The bound- 
aries of the region in the parameter space permitting transition of number of critical 
points in this case are associated with saddle-centre bifurcation or merging of a pair 
of roots of the equation (EqfTTll. Now all these equations are polynomial equations. 
As a general rule the discriminant of a polynomial, 

P„{x) = a„x" + a„_ix""^ H h aix + oq, (24) 

can be expressed as in terms of its roots, Xi S, SLS 

The discriminant may be expressed as the determinant of a matrix called the Sylvester 
matrix (see, e.g., |http://mathworld.wolfram.com/PolynoriuaLDiscriminant.html| and 
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references therein), 



a„ a„_i a„_2 ••• ai oo 0... 

^ _ ... a„ fln-i • ■ ■ oi 10 

7ia„ (71 — l)a„_i (rt — 2)a„_2 ••• loi ... ... 

nUn (n — l)a„_i (tt, — 2)a„_2 ••• loi ... 



... nan (n — l)an-i (n — 2)an-2 ■■■ loi 

(26) 

up to a factor 

Putting n — 8, det (S) will be zero on the above mentioned boundaries and 
actually it is so. Here the plot of det {S) = for the polytropic flow (i.e. for the poly- 
nomial in Tc in EqfTTli in £-X space is shown in Fig|2] The curve exactly conforms 
with the corresponding boundary curve in FigJT] drawn on the basis of the previous 
method. So this procedure may be thought of as a much easier alternative to find the 
multi-critical parametric values; though this method cannot give the exact number of 
critical points in each region of the parameter space. 



5 Discussion 

Our methodology is based on the algebraic form of the first integral obtained by solv- 
ing the radial momentum equation (the Euler equation to be more specific, since we 
are confined to the inviscid flow only). The structure for such a first integral has to be 
a formal polynomial with appropriate constant co-efficients. For general relativistic 
accretion in the Kerr metric, the expression for the energy first integral can not be re- 
duced to such a polynomial form (see, e.g., |33 1 for the detail form of such algebraic 
expression). Hence, the Sturm's generalized chain can not be constructed for such 
accretion flow. Alternative methodology are required to investigate the multi-critical 
behaviour for such kind of accretion. 

Using the method illustrated in this work, it is possible to find out how many 
critical points a transonic black hole accretion flow can have. It is thus possible to 
predict whether such accretion flow can have multi-critical properties for a certain 
specific value/domain of the initial boundary conditions. It is, however, not possible 
to investigate, using the eigenvalue analysis as illustrated in [8','T2'|, the nature of 
such critical points - i.e., whether they are of saddle type or are of centre type, since 
such prediction requires the exact location of the critical points (the value of the 
roots of the polynomial). However, the theory of dynamical systems ensures that no 
two consecutive critical points be of same nature (both saddle or both centre). On 
the other hand, our experience predicts (it is rather a documented fact) that for all 
kind of multi-critical black hole accretion, irrespective of the equation of state, the 
space time geometry or the flow configuration used, one has two saddle type critical 
points and one centre type critical point flanked by them (see, e.g., |[33l and ll34l for 
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farther detail). Hence if the application of the Sturm's generalized chain ensures the 
presence of three critical points, we can say that out of those three critical points, 
accretion flow will have two saddle type critical points, hence a specific subset of the 
solution having three roots corresponding to the first integral polynomial, can make 
transonic transition for more than one times, if appropriate conditions for connecting 
the flow through the outer critical point and for flow through the inner critical points 
are available, see, e.g., 133 1 for further discussion. 

In this work we have considered only inviscid accretion. Our methodology of in- 
vestigating the multi-critical properties, however, is expected to be equally valid for 
the viscous accretion disc as well. For the viscous flow, the radial momentum con- 
servation equation involving the first order space derivative of the dynamical flow 
velocity will certainly provide a first integral of motion upon integration. Because 
of the fact that a viscous accretion disc is not a non-dissipative system, such con- 
stant of motion, however, can never be identified with the specific energy of the flow. 
The integral solution of the radial momentum equation would then be an algebraic 
expression of various flow variables and would perhaps involve certain initial bound- 
ary conditions as well. Such an algebraic expression would actually be a constant of 
motion. What exactly would that expression physically signify, would definitely be 
hard to realize. However, one may perhaps arbitrarily parameterize that conserved 
algebraic expression using some astrophysically relevant outer boundary conditions, 
and if such algebraic expressions can finally be reduced, using the appropriate critical 
point conditions, to an algebraic polynomial form of the critical points, construction 
of a generalized Sturm chain can be made possible to find out how many critical 
points such an accretion flow can have subjected to the specific initial boundary con- 
dition. Since for accretion onto astrophysical black holes, having multiple critical 
points is a necessary (but not sufficient) condition to undergo shock transition, one 
can thus analytically predict, at least to some extent, which particular class of viscous 
accretion disc are susceptible for shock formation phenomena. 

Our work, as we believe, can have a broader perspective as well, in the field of 
the study of dynamical systems in general. For a first order autonomous dynamical 
system, provided one can evaluate the critical point conditions, the corresponding 
generalized nth degree algebraic equation involving the position co-ordinate and one 
(or more) first integral of motion can be constructed. If such algebraic equation can 
finally be reduced to a nth degree polynomial with well defined domain for the con- 
stant CO efficient, one can easily find out the maximal number of fixed points of such 
dynamical systems. 
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6 Appendix - 1 : Proof of the Sylvester's theorem: 

First note that the Sturm sequence (/o, ■■■fk) is (up to signs) equal to the sequence 
obtained from the EucHdean algorithm. Define a new sequence (go, 9k) by Qi = 
fi/fk for i & {0,...,k}. Note that the number of sign changes in {fo{x), fi{x)) 
(resp. fi{x), fi+i{x))) and the number of sign changes in {go{x),gi{x)) 

(resp. {gi^i{x), gi{x), gi^i{x))) coincide for any x which is not a root of /. Note 
also that the roots of go are exactly the roots of / which are not roots of g. Observe 
that fori e 0, ...,k,gi-\ and^fj are relatively prime. We consider, now, how v{f,g;x) 
behaves when x passes through a root c of a polynomial g^. If c is a root of go, then 
it is not a root of gi. We write /'(c) > (resp. < 0) if /' is positive ( resp. negative 
) immediately to the left of c. The sign of /'(c+) is defined similarly. Now we recall 
the following result: if i? is a real closed field, / G a,b E R with a < b and 

if the derivative /' is positive (resp. negative) on ]a, b[, then / is strictly increasing 
(resp. strictly decreasing) on [a, b]. Then, according to the signs of g{c), f'{c-) and 
f'{c+) we have the following 8 cases: 

5(c)>0,/'(c_)>0,/'(c+)>0 





C- 


c 




/ 







+ 


fg 


+ 




+ 



fl(c)<0,f(c_)>0,/'(c+)>0 







c 




/ 







+ 


fg 









5(c) >0,/'(c_) <0,/'(c+) >0 







c 




./ 


+ 





+ 


I'g 






+ 



5(c)<0,/'(c_)<0,/'(c+)>0 







c 




/ 


+ 





+ 


fg 


+ 







fl(c)>0,/'(c_)>0,/'(c+)<0 
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C- 
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c+ 


/ 









fg 


+ 







5(c)<0,/'(c_)>0,/'(c+)<0 







c 




/ 









fg 






+ 



5(c)>0,/'(c_)<0,/'(c+)<0 







c 




/ 


+ 







/'5 









5(c)<0,/'(c_)<0,/'(c+)<0 







c 




/ 


+ 







fg 


+ 




+ 



In every as x passes through c, the number of sign changes in {fo{x), fi{x)) 
decreases by 1 if g(c) > 0, and increases by 1 if g(c) < 0. If c is a root of with 
i = 1, .../c, thenitis neitherarootof (7i_i norarootof gi+i, and (?i_i(c)(?j+i(c) < 0, 
by the definition of the sequence. Passing through c does not lead to any modification 
of the number of sign changes in (a;), fi{x), /j+i {x)) in this case. 

Proof of the Sturm's theorem: Using = 1 in previous theorem. 
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7 Appendix - II: Explicit expressions for the co-efficients for the Sturm chain 
constructed for the relativistic axisymmetric accretion 
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62^5 , 63^5 

73 = + ( —je3 - ds 
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61^5 d4 63^5 
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